A Real-Space Full Multigrid study of the fragmentation of Li^ clusters. 
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We have studied the fragmentation of Li^ clusters into the two experimentally observed products 
Lijj" + Li2 and Li^~ + Li. The ground state structures for the two fragmentation channels are found 
by Molecular Dynamics Simulated Annealing in the framework of Local Density Functional theory. 
Energetics considerations suggest that the fragmentation process is dominated by non-equilibrium 
processes. We use a real-space approach to solve the Kohn-Sham problem, where the Laplacian 
operator is discretized according to the Mehrstellen scheme, and take advantage of a Full MultiGrid 
(FMG) strategy to accelerate convergence. When applied to isolated clusters, we find our FMG 
method to be more efficient than state-of-the-art plane wave calculations. 



PACS numbers:66.10Ed,71.55Jv,72.20Jv 
I. INTRODUCTION 

In recent years a number of real-space grid-based algo- 
rithms have been proposed to study the electronic prop- 
erties of condensed matter systems. Most of these ap- 
proaches, based on Density Functional (DF) theory, eval- 
uate in real space the action of the effective potential and 
of the kinetic energy operators that enter the one-electron 
DF equations. 

When based on an accurate approximation of the 
Laplacian operator, real-space methods exhibit sev- 
eral advantages when compared to Plane Waves (PW) 
schemes. In particular: 

a) Boundary conditions are not constrained to be pe- 
riodic, thus facilitating the study of charged systems or 
of systems having finite electric multipole moments (e.g. 
isolated clusters). 

b) In real space, it is technically easier to increase the 
resolution of the grid locally. This feature is particularly 
appealing in the case of first-row and transition-metal el- 
ements and for the study of clusters: a higher density 
of grid points is usually needed in the regions where the 
ions are located, whereas a lower density is sufficient on 
the vacuum side to describe the exponentially small tails 
of the electronic wavefunctions. The use of locally re- 
fined meshes may lead to a significant reduction of the 
computational effort. 

c) With a suitable discretization of the Laplacian op- 
erator, all the operations needed to compute the total 
energy of the system are short ranged, hence facilitating 
the parallelization of the computer code. 

Among the various approaches reported in the litera- 
ture, let us mention the use of high-order finite differ- 
ence methods for representing the laplacian, combined 
with the use of soft non-local pseudopotentials on uni- 
form grids jl| to calculate the electronic structure and the 
short time dynamics of small molecular systems. Gygi 
et al. extended the real-space adaptive-coordinate 
method of Ref. || to perform LDA electronic structure 
calculations of small molecules. Ref jlj reports a similar 
approach. Other recent real-space methods include the 
use of wavelet basis sets 15101 , the finite- elements method 



|@,|| and the Multigrid (MG) calculations of first and 
second row atoms Jsjl . 

Briggs et al. have recently proposed an efficient 

method that combines an accurate discretization of the 
Laplacian (Mehrstellen) with a MG acceleration scheme 
to solve Poisson and Kohn-Sham (KS) equations in large- 
scale LDA calculations. 

Our real-space calculations follows Briggs' et al. [po|JTT] 
in the use of the (Mehrstellen) discretization scheme, but 
differs from previous calculations in a significant way in 
the scheme chosen to accelerate convergence. Whereas in 
Ref. jllj a MG solver was used to solve Poisson and KS 
equations we implement here a Full MultiGrid (FMG) 
diagonalization procedure. 

Given some similarities between our method and the 
one of Ref. JO]] , we will concentrate on their differences 
and refer the interested reader to Ref. [^o|Jll| l for technical 
details common to both approaches. 

This paper is organized as follows: in Section II we dis- 
cuss our Full Multigrid (FMG) algorithm, in Section III 
we demonstrate its accuracy and performances by run- 
ning a few tests on simple molecules. We then apply our 
scheme to the study of the fragmentation of Li^ clusters 
in Section IV. Section V contains a summary and some 
concluding remarks. 



II. COMPUTATIONAL METHOD 

A. Multigrid Methods 

MG methods were introduced in the 70 's by Brandt 
Jl2[ as a tool to solve elliptic PDE discretized on a grid 
of TV points in O(N) operations. His idea was to use 
some auxiliary set of grids to speed up the convergence 
of straight relaxation methods [13] . In these methods the 
value at the mesh points are repeatedly updated accord- 
ing to the discretized differential equation. This process 
enforces a local consistency between the updated value 
and that of its neighbors. By performing enough re- 
laxation cycles the exact solution propagates from the 
boundaries, where it is fixed by the boundary conditions, 
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to the interior of the grid. Relaxation iterations reduce 
quickly the components of the error with wavelength 
comparable to the grid spacing but remain quite inef- 
fective in reducing error components with larger wave- 
length. The MG strategy aims at removing these low 
frequency errors by noticing that a) low frequencies on 
a fine grid become higher frequencies on a coarser grid, 
and b) the error on the solution of a PDE with fixed 
boundary conditions satisfies a similar PDE, with zero 
boundary conditions. Solving this PDE for the error on 
a coarse grid allows to remove the low frequency errors 
that spoil the fine grid solution. This two-level strategy, 
often referred to as "V-cycle" can naturally be 

applied recursively and involve an arbitrary number of 
levels. These "V-cycles" allow to converge very rapidly 
to the right solution. They enhance drastically the effi- 
ciency of relaxation methods and are referred to as the 
MG method. The interested reader is referred to Ref. 
|jl4| and to references therein for more details. 

MG methods are particularly well adapted to solve 
Poisson equation and to large scale eigenvalue problems 

HQ- 

An improvement over MG methods is made possible by 
using a FMG strategy. While both approaches make use 
of Brandt's idea and use the auxiliary set of grids to per- 
form "V-cycles", the FMG approach is somewhat more 
complete in that it also uses these grids to construct a 
good first guess at almost no cost. In other words, while 
MG codes start the computation on the final grid and re- 
sort to "V-cycles" to speed up convergence, FMG codes 
solve the problem on some coarse grid first and, once 
the solution is converged well enough, expand it on the 
next finer grid. The coarse grid solution is a precondi- 
tioning for the initial guess on the next finer grid. This 
preconditioner is used up to the finest grid. FMG codes 
take obviously advantage of the "V-cycles" on each of the 
intermediate levels. 

As a consequence, FMG schemes further enhance the 
rate of convergence of MG algorithms since they reduce 
the residuals at each iteration In practice, FMG is 
only slightly more complex to implement than MG. 



B. Grid representation of the Kohn-Sham problem. 

In the Kohn-Sham [[l7| formulation of Density Func- 
tional Theory the total energy E[{^i}, {R a }] of a 
system of N a ion cores located at {R a } and of N elec- 
trons writes (in atomic units and considering only doubly 
occupied states): 

OCC r. 

£[{^},{R a }] = 2]T / #[-~V 2 ]Vidr 

1=1 

p(r)p(r>) 



E xc [p(r)} + E, 
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where the {ipi}'s are the Kohn-Sham (KS) electronic or- 
bitals. 

The first, second and third terms of the right hand 
side are the kinetic, the electron-electron (Hartree) and 
the electron-ion energies, respectively. The fourth term 
stands for the exchange-correlation contribution, while 
the last one corresponds to the ion-ion electrostatic re- 
pulsion. The electron density is given in terms of the 
N/2 occupied KS orbitals by: 



OCC 

p(r) = 2j2Mr)\ 



where the KS orbitals are the solutions of: 
[-iv 2 + V e //]^(r) = e i ^(r), 



(2) 



(3) 



subject to the orthonormality constraint < ipi\ipj >= 
The effective potential V e ff (r) is given by: 



V eff (r) = V lon (r) + Vff(r) + V xc (v) 

P(r') 



(4) 



where Vff(r) = J \ ^._J \ dvi is the Hartree potential and 
V xc (y) = SE xc /5p(r) the exchange-correlation potential. 
In our calculations the electron-ion interaction Vi on is 
represented by ab initio norm-conserving pseudopoten- 
tials. We use the Kleinman-By lander |l9| form in real 
space |20| for the non-local part of Vi on . 

For simplicity we reported the equations for param- 
agnetic systems i.e. for systems with all states doubly 
occupied. The generalization of Eqs. (|l|) and (||) to spin- 
polarized system (Local Spin Density, LSD) is straight- 
forward |2lf| . For our LDA and LSD computations we 
resort to the Perdew and Zunger [ p2[ parametrization of 
the Monte Carlo data of Ceperley and Alder 23 1 . 

Following Ref. |Ioj , we discretized the Laplacian opera- 
tor in Eq. (3) according to a Mehrstellen scheme of order 
0(h 4 ) |24|. The latter was indeed proven to provide at 
least the accuracy of schemes using a sixth-order central 
finite-difference discretization of the Laplacian |]Io| , |rT| , 
with the advantage of using more local data (i.e. only 
up to second-nearest neighbors for a cubic mesh). 

The use of a uniformly spaced grid with spacing h al- 
lows to define an effective energy cut-off equal to that of 
the PW calculations that use the same real-space grid for 
their FFT's. Assuming that the electronic charge density 
in the plane wave code is expanded with an energy cut- 
off four times as large as that of the wavefunctions, we 
define G 2 max = tt 2 /4/i 2 Ry. 

In our code all integrations are performed according to 
the three-dimensional trapezoidal rule: 



J g(r)dr = h 3 J29«i,j,k)) 

ijk 



(5) 



(r)V ion (v)p(r)dr 



We notice that for high accuracy, it is essential that all 
the integrands g(r) be band- limited in the sense that 
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their Fourier transform must have zero magnitude in the 
frequency range G > G max = n/2h. While this condi- 
tion is automatically fulfilled in PW calculations, since 
the basis set is cut-off at a specific plane- wave energy, in 
real-space calculations high-frequency components above 
the natural cut-off G\^ ax can manifest themselves on the 
grid. In particular, if the pseudopotential Vi on contains 
high frequency components near or above G max , these 
high frequency components may be aliased to lower fre- 
quencies in an unpredictable manner | jTo[ leading to un- 
physical variations in the total energy as the ions move 
with respect to the grid. 

This defect can be overcome by explicitly eliminating 
the high frequency components of the pseudopotentials 
by Fourier filtering, as shown in the context of PW by 
King et al. p5| ] . In this work, we followed their prescrip- 
tion and removed from the non-local pseudopotential the 
Fourier components near G m ax- This filtering procedure 
doesn't introduce any significant computational overload 
since it is done once for all at the beginning of the calcu- 
lation. In real space every update of the non-local part 
of the total energy and of ionic forces scales as the square 
of the number of atoms in the unit cell, as opposed to the 
scaling as the cube of the system size in the conventional 
reciprocal-space formulation p5| . 

C. Full MultiGrid approach to the solution of the 
Kohn-Sham problem. 

In order to determine the ground-state electronic den- 
sity for a given configuration of ionic cores, we solve self- 
consistently Eq. (||) using a FMG approach along the 
lines described below. 

In our calculations we use a cubic simulation cell of 
volume L 3 with k max uniform samplings. Calling k — 
0, 1, ...kmax the level number, the grid at level k contains 
N% points spaced by h k = L/(N k -l), with N k = 2 k (n - 
1) + 1. In most cases we use three grid levels k = 0, 1,2, 
refered to as "coarse", "intermediate" and "fine" in the 
following. On the "coarse" and "intermediate" grids we 
take Vi on as purely local and approximate it with the s- 
componcnt of the non-local pseudopotential. The choice 
of the coarsest grid and thus of no, is dictated by the 
constraint that it must contain enough points to allow for 
a meaningful representation of all the KS orbitals we are 
interested in: i.e. it must allow for the orthogonalization 
of the KS orbitals and still leave them enough degrees of 
freedom to build up a reasonable density. 

On every grid, the relaxations are of the Gauss-Seidel 
type flij and are performed according to the red-black 
ordering. The restriction and interpolation procedures 
used in the MG-cycles are the projection of a weighted 
average and a simple tri-linear interpolation. The averag- 
ing procedure involves the first, second and third nearest 
neighbors of the fine point to be transferred. 

Our FMG solver starts on the "coarse" level of Uq 
points, where we perform a few self-consistent iterations 



only. The KS orbitals are then transferred to the "in- 
termediate" grid. Once the orbitals are converged well 
enough they are interpolated on the "fine" grid. We solve 
on the latter the final KS problem. On the intermediate 
and fine grid we take full advantage of the MG strategy 
we will describe in the following. 

We may identify here a first fundamental difference be- 
tween our approach and that of Briggs et al. : while 
they solve the KS problem directly on the fine grid and 
resort to the auxiliary set of coarse grids in a MG strat- 
egy, we, in contrast, implement a FMG solver in order to 
precondition the initial KS orbitals on the fine grid and 
thus start our computations on the coarsest grid. Possi- 
ble advantages of our approach are discussed in Section 
III. We will outline another important difference in sub- 
section 2 below. 



1. Poisson equation 

The Hartree potential Vh in Eq. (4) is the solution of 
the Poisson equation: 

\7 2 V H = -47rp (6) 

where p is the electronic charge density. 

Since we are interested here in isolated clusters, we 
determine the boundary conditions on the potential Vh 
from the multipoles of the charge density. When the sim- 
ulation cell is large enough, this is a good approximation 
that allows the study of charged as well as of neutral sys- 
tems. Other boundary conditions (e.g. periodic) can be 
imposed without difficulty Jl(j] . 

Our FMG Poisson solver does usually have more levels 
then the KS-solver. The coarsest grid for Poisson has 
n 0P points and nop relates to the no of the KS problem 
by no = 2 J (nop — 1) + 1 with j such that nop is the 
smallest possible integer larger than one. On the coarsest 
grid we use an exact solver, a direct matrix inversion if 
necessary. The determination of Vh on the set of coarse 
grids is of course possible since we know p on a fine grid. 
The "V-cycles" for Poisson are based on the following set 
of equations Q : 

V H = V H + e (7) 
r = 4np + V 2 V H (8) 

V 2 e = r (9) 

e stands for the error on the actual solution Vh and the 
residual r is a measure of this error. We can determine 
e by solving Eq.(9) for zero boundary conditions. FMG 
approaches allow to converge the solution of the Pois- 
son equation to the desired accuracy with a minimum 
number of V-cycles on the finest grid. Notice that the 
addition of the extra coarse grid levels is computation- 
ally inexpensive since at each level the number of grid 
points is reduced by a factor of approximately eight. 
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2. Kohn-Sham equations 

Solving the KS equations (Q) is obviously far more com- 
plex than solving the Poisson problem since Eqs. (||) are 
non-linear: wavefunctions and eigenvalues must be de- 
termined simultaneously. In the spirit of the FMG pro- 
cedure, we start our computations with some guess on 
the "coarse" grid, and get there a solution as accurate as 
possible. 

Since we already explained the basics of our FMG ap- 
proach to the KS problem, we just state here the equa- 
tions our "V-cycles" are based on. In the following ipi 
stands for the actual solution, ipi for the exact solution, 
and r for the residual associated to ipi. 



ipi = aipi + ip i± 



-i>i + (v- a)^ 



(10) 



(11) 



(12) 



We determine an approximate ipij_ by a few relaxions 
according to Eq. jl2] ) and use Eq.(^p|) to improve ipi. 

We evaluate the eigenvalue associated to {ipi} by re- 
sorting to Rayleigh's formula |l^| and ensure the orthog- 
onality of the states with a Gram-Schmidt procedure. 

In order to avoid instabilities due to charge "sloshing" , 



we mix the new and old densities: p r , 



(1 - a)p id + 



2aJ2i=i W: where < a < 1 (usually a value of a be- 
tween 0.3 and 0.5 is used in our calculations). 

At this point we stress a second significant difference 
between our scheme and that of Briggs ct al. 11 1 . While 
we base our V-cycle on (12) they resort to the Poisson 
equation: 



(13) 



implicitly assuming that the leading terms entering the 
residual r are due to an overestimate of the kinetic energy 
contribution. Furthermore, while we use the Mehrstellen 
discretization on every grid level, they use it on their 
finest grid only. On the other grids they approximate 
the Laplacian operator with a 7-point central finite- 
difference. 

To conclude this section, we note that, as suggested by 
Briggs et al. J2^jll[], it is crucial to perform a Ritz pro- 
jection every now and then during the self-consistent 
iterations in order to unmix eigenstates that may be close 
in energy. This procedure improves the convergence to 
the exact self-consistent solution. 



III. TESTS AND RESULTS 

In order to check the accuracy and robustness of our 
approach and to compare its efficiency with conventional 



PW calculations we made a number of tests on small 
molecules. 

First of all, we verify that the spurious effects due 
the components of the pseudopotential with frequencies 
higher than G max = n/2h are actually made very small 
when the pseudopotentials are filtered according to the 
procedure p5] | mentioned in Section 2.B. These spurious 
effects usually appear as an unphysical dependence of the 
total energy on the positions of the ions relative to the 
grid points. 

To check this we consider an isolated C*2 molecule in 
the middle of a cubic cell and assume zero boundary con- 
ditions on the wavefunctions. We make a series of cal- 
culations where the C2 dimcr is displaced rigidly relative 
to the grid along a direction parallel to the cell edge and 
evaluate the total energy for each position. The largest 
displacement from the center of the cell is h/2, where 
h is the grid spacing of the finest grid used in the cal- 
culation. We assume a cubic cell of side L — 14a.u., 
which we checked to be large enough to justify our zero 
boundary conditions, and use three uniform grids of 17 3 , 
33 3 and 65 3 points. The Poisson equation is solved us- 
ing additional coarser grids down to a mesh of 3 3 points. 
The finest grid used for the KS problem corresponds to a 
plane-wave energy cut-off E cu t = 52 Ry, which is suf- 
ficient to give reasonably well converged properties of 
Carbon systems. We use ab-initio norm-conserving pseu- 
dopotentials J27|] with s-non locality. 

Usually about 20 self-consistency cycles are sufficient 
to achieve a good convergence for a given dimer posi- 
tion. Every of these cycles consists of 2 V-cycles, all 
the relaxations are repeated twice and the eigenvalues 
are updated after every V-cycle. We perform the Ritz 
projection, mentioned in the previous Section, every 5 
self-consistent cycle. 

We find that the variation of the total energy as the 
dimer is displaced across the cell is at most 0.5 meV, 
showing that aliasing effects are indeed very small, at 
least on the scale of energies we are interested in. 

We have also calculated the equilibrium distances 
and the vibrational frequencies of a number of simple 
molecules (C2, O2, CO, Sia) and compared our results 
with experiments. When possible, comparisons were also 
made with the results obtained with a PW code using an 
energy cut-off corresponding to the mesh used in our real- 
space calculation. The two sets of results are in very good 
agreement with each other since we didn't find discrep- 
ancies exceeding 1% both for the equilibrium distances 
and for vibrational frequencies. 

As an example, we report the results of calculations for 
a C2 dimer (those for other molecules are very similar). 
The upper pannel of fig. (|l|) displays the total energy Eq. 
(0) of the C 2 dimer at 5 different bond length (squares) 
along with a fourth order polynomial fit (solid line). 
From the fit we find the equilibrium distance do = 1.26 A 
and the vibrational frequency u) v = 1860 cm -1 to be 
compared with the experimental values do = 1.24 A and 
io v = 1854 cm -1 . In the lower panel of Fig. (nh, we 
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compare the forces, calculated according to Hellman- 
Feynman's theorem (squares) to the analytical deriva- 
tive (solid line) of the polynomial curve that interpolates 
the total energy values. The quality of the force compu- 
tation can be judged from the perfect matching of the 
points with the solid line. 
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FIG. 1. Upper panel: Total energy of the Ci dimer as 
a function of the C — C separation. The squares are the 
calculated points, the solid line is a fit with a 4 — th order 
polynomial. Lower panel: Force acting on the C atoms. The 
squares are the calculated values, the solid line is obtained by 
derivative of the analytic curve shown with a solid line in the 
upper panel. 

For the systems we considered, we found the speed of 
convergence of the FMG to the electronic ground-state, 
for a given ionic configuration, substantially increased 
with respect to that of state-of-the-art Car-Parrinello 
(CP) like methods 28 1, i.e. much shorter CPU times 
are required to achieve the same accuracy in the self- 
consistent charge density, total energy and ionic forces. 
For the purpose of comparison we ran a CP code based 
on a Damped Molecular Dynamics relaxation of the 
electronic degrees of freedom, with the maximum times 
step allowed to ensure stability during the minimization. 
More efficient CP schemes exist p9| , where the conver- 
gence of the total energy and forces is accelerated by 
preconditioning techniques, which allow larger integra- 
tion steps than those usually required in standard CP 
codes. Even in the latter case our FMG code shows bet- 
ter performances in minimizing the electronic energy at 
fixed atomic positions, at least for isolated molecules or 
clusters. Furthermore, as shown in the following, our ap- 
proach allows to use larger time steps in Molecular Dy- 
namics, and thus it is preferable for Molecular Dynamics 
simulations. 
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FIG. 2. Comparison between the convergence rates of 
the total energy and forces for the C2 molecule for our FMG 
calculations (lower solid and dotted curves respectively) and 
for plane- wave calculations (upper solid and dotted curves). 

An example of the convergence rate of the total en- 
ergy and forces is reported in Fig. (^J) for the case of a 
C2 dimer. The logarithm of the deviation of the total 
energy E and force F from their values at convergence, 
E and Fn, are shown as a function of the number of self- 
consistency cycles (lower horizontal axis) for our FMG 
(lower two curves) and as a function of the number of 
time steps (upper horizontal axis) for the CP run (upper 
two curves). Both computations required similar CPU 
time. 

We did not attempt a direct comparison between our 
FMG and the MG approach chosen by Briggs et al. [|l0| , 
although we expect similar performances as far as the re- 
laxation rate is concerned. We note however that, since 
our FMG calculations are started on the coarsest grid, 
the use of a very rough initial guess for the electronic 
wavefunctions (random numbers) and density (superpo- 
sition of atomic densities) is usually sufficient to start 
the computation. More educated guesses do not help in 
reducing the total CPU time required to achieve conver- 
gence. On the other hand, when we use a MG schedule, 
i.e. start on the finest grid, a very good initial guess 
for the electron wavefunctions is mandatory to achieve 
self-consistency with the same CPU time. 

Having established the accuracy of our method for the 
calculation of energy and forces, we checked its perfor- 
mances in the search of the lowest energy structure of 
a small metallic cluster, AIq. We use a nonlocal pseu- 
dopotential, with s-non locality. The side of our cell is 
L = 25 a. u.. It is sampled with three different grids of 
13 3 , 25 3 and 49 3 points. The effective energy cut-off on 
the latter grid is 9 Ry, which is sufficient to represent the 
pseudo charge density of Al atoms in a condensed phase 
pOfl. Our starting geometry is arbitrary, i.e. a planar 
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hexagonal ring of Al atoms. To search for the lowest 
energy structure of this cluster, we used a Simulated An- 
nealing schedule where the temperature is initially raised 
to T ~ 1000 °K by rescaling the atomic velocities and 
then slowly reduced to T = 0. The history of the MD 
run is recorded in Fig. (^|). The cluster temperature is 
shown as a function of the number of MD time steps. 
The time step is At = 200 a.u. and the total simula- 
tion time is ~ 3ps. After a first annealing cycle (first 
~ 200 time steps in Fig. (||)) the cluster got caught 
in a local minimum of the total energy, corresponding 
to an open, quasi two-dimensional, distorted structure, 
~ 0.17 eV I atom higher in energy than the stable struc- 
ture. We thus reraised the temperature, annealed, and 
performed a Steepest Descent relaxation when in sight 
of the ground-state structure. The resulting geometry is 
displayed in Fig. (||). This structure is characterized by 
two bond lengths, l\ — 2.58 A and li = 2.97 A, and ex- 
hibits Z?3rf symmetry, in agreement with previous LDA 
calculations p(J. 



A 




200 400 600 

§ of time steps 



FIG. 3. Simulated annealing cycle for the Ale cluster. In 
the last 30 steps of the run a Steepest Descent relaxation on 
the ionic coordinates of the cluster has been performed. The 
final structure is shown in the inset. Bonds AE, AF, BE, 
CF, BD and CD in the lowest energy structure have length 
h = 2.58 A, while bonds AB, AC, BC, EE, DE and DF 
have length h = 2.98 A. 

IV. THE FRAGMENTATION OF LI+ CLUSTERS 

After checking the accuracy and robustness of our im- 
plementation of the FMG scheme, we have applied it to 
the problem of Li^ fragmentation, which has been the 
subject of recent experiments |3l[| . 

Experimental studies of fragmentation of small metal- 
lic clusters provide useful information for understanding 



the properties of matter at the small aggregation limit 
and the size-dependent evolution towards bulk behavior. 
In particular, the dynamics of unimolecular clusters dis- 
sociation characterizes the nature of the energy partition- 
ing among the internal vibration modes and its under- 
standing may provide important information on the clus- 
ters properties. Statistical models of energy partitioning 
allow, for instance, the determination of the binding en- 
ergy of the clusters from their dissociation rates. Cluster 
dissociation processes may also give some insight into the 
important problem of estimating the "melting" tempera- 
ture of the cluster, i.e. the temperature above which the 
cluster becomes "liquid" i.e. with no well defined struc- 
ture. Experimentally, the temperature dependence of the 
heat capacity has been measured, for instance, from the 
photofragmentation mass spectrum of alkali clusters al- 
lowing the identification of the cluster melting point |$2| . 

In most cases, neutral and singly charged metal clus- 
ters are intrinsically stable and their dissociation is en- 
dothcrmic. The excess energy necessary to promote dis- 
sociation is provided by the photon energy in photoab- 
sorption experiments. 

Evaporation of ionized lithium clusters has been stud- 
ied recently in experiments where high-power lasers were 
used to ionize and excite clusters, and where the evap- 
oration into smaller products took place in a time long 
with respect to the timescale set by the vibrational fre- 
quencies of the cluster. A competition between the evap- 
oration of monomers and dimers, with the production of 
Li^ and Lig cations respectively, has been observed |n] . 
The dynamic behavior was found to depend critically on 
the ionization conditions. The evaporation of dimers was 
found to be the most significant fragmentation channel 
under irradiation with a AeV laser, i.e. slightly above the 
photoionization threshold, whereas suddenly both evap- 
orative paths appeared as soon as the energy exceeded 
<l.3eV and found to be of comparable intensity when the 
cluster was ionized with a A.5eV laser. 

If one assumes that during the photodissociation the 
internal energy due to multi-photon absorption is ran- 
domly distributed in the metastable cluster over the 
3n — 6 internal modes, then the dissociation of a Li+ 
cluster occurs as soon as enough internal energy becomes 
localized in a single mode, so as to overcome the fragment 
binding energy D^: 

D+(l) = £7(ii+_ 1 ) + E(Li) - E(Li+) (14) 

or 

D+ (2) = E(Li+_ 2 ) + E(Li 2 ) - E(Li+) (15) 

for the two possible fragmentation channels, respectively. 
Higher dissociation rates are expected for the channel 
characterized by a lower dissociation energy. From the 
analysis of the experimental results pl[ , however, non- 
statistical effects seem to play a major role in the pro- 
cess of fragmentation of Li^, i.e., even a long time after 
the photoexcitation process, the internal energy appears 
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not to be randomly distributed among all the vibrational 
modes (similar non-statistical effects have been observed 
for Na^ clusters [[33|). As a consequence, dissociation 
may occur along less energetically favorable channels. In 
fact, as we will show in the following, our total energy 
calculations predict slightly lower dissociation energy for 
the dissociation channel Lif ± — > Lii + Li, with the pro- 
duction of monomers. 

Configuration Interaction (CI) calculations of the ge- 
ometry and electronic structure of small neutral Li n and 
cationic Li+ clusters, with n up to 10, have been reported 
in the literature Q]. The binding energies £)+(l), char- 
acterizing the dissociation channel Lz+ — > Li^ l _ 1 + Li, 
exhibit alternating maxima and minima for clusters with 
odd or even number of atoms respectively, showing their 
smaller or larger stability for this dissociation channel. 
For instance, smaller stability characterizes the dissocia- 
tion of Li\ , Lig ,Lif with the production of monomers. 
A less characteristic behavior is found in D+(2) for the 
Li+ — » Li^_ 2 + Lii dissociation process. 

Whereas, on the basis of CI calculations |$4j, small 
clusters (with 3 to 6 Li atoms) can be considered as de- 
formed sections of the (111) plane in the fee crystals, the 
structure of larger clusters can be described (although 
with a few exceptions) as more open structures, com- 
posed of deformed tetrahedrons appropriately sharing 
their triangular sides. 

The experimental binding energies for Li n cations were 
derived from the interpretation of photodissociation ex- 
periments |p5| . They are determined from a statistical 
treatment of unimolecular decay rates, where the energy 
is assumed to be statistically distributed over the various 
nuclear degrees of freedom. In Rcf. |35| the binding ener- 
gies of the cationic Li+ clusters with n up to 10 obtained 
by this phenomenological method are compared with re- 
sults of CI calculations. The discrepancy between CI and 
statistical cluster dissociation energies can however be as 
large as 0.4 eV. 

The same statistical model predicts |3j| the values 
D+(l) = 1.30 and D±(2) = 1.07 eV for the two 
evaporative channels of Li{ x clusters, but no first prin- 
ciples calculations for this cluster exist to support these 
estimates. We provide here results, obtained with our 
FMG method, for the dissociation energies of Li^ x clus- 
ters along their two main fragmentation paths. At vari- 
ance with the results |35|] quoted above, we find that the 
fragmentation into monomers is slightly favored on ener- 
getics grounds, i.e. Z?n(l) + is lower than £>n(2) + . 

The Full Multigrid scheme is implemented on a set of 
three grids of N 3 points, with N = 21,41,81. We used 
the LSD version of Eq. (Q), as described in Section 2. 
The electron-ion interaction is modeled with a nonlocal 
pseudopotential |27]] , with s-non locality. The side of the 
cell is L = 36a.u.. This means that on the finest grid 
the effective energy cut-off is 12 Ry, which is sufficient 
to represent the pseudo charge density of Li atoms in a 
condensed phase. The values d eq — 2.65 A (2.67) and 
uj v = 348 cm -1 (351), obtained for the equilibrium dis- 



tance and the vibrational frequency (in parenthesis we 
report the experimental values) of the Li 2 dimer show 
the reliability of the pseudopotential to provide an accu- 
rate description of the system. 

An unbiased search, Simulated Annealing, is used to 
obtain the lowest energy atomic configurations for the 
cations Li+, with n = 9,10 and 11. Annealing is typically 
started from a randomly generated initial configuration 
for the ionic cores; the "temperature" (defined in terms 
of the total ionic kinetic energy) is initially raised to T ~ 
800 — 1000 °K by rescaling the atomic velocities and then 
slowly reduced to T = 0. We used a simulation time step 
At = 200 a.u.. 

The resulting structure for Li\ x , very similar to the 
lowest energy structure of the neutral Lin cluster, is 
shown in Fig. (^). The cluster has C^v symmetry and the 
relevant bond lengths are reported in the Figure (for com- 
parison, we remind that the experimental interatomic 
distance in the fee Li crystal is 3.10 A). 




FIG. 4. Equilibrium geometry of the Li{ x cation. Bond 
lengths are shown in A. 

We also calculate the dipole polarizability of the clus- 
ter by applying a small uniform electric field E and by 
computing the resulting dipole moment P. The field in- 
tensity must be small enough to produce a dipole moment 
proportional to the field. The polarizability is then eval- 
uated asa = P/E. We find at\\ = 128 A 3 and a± = 76 A 3 
for the polarizabilities along the cluster axis connecting 
the two cap atoms at the far ends of the molecule, and 
perpendicular to it. 
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FIG. 5. Equilibrium geometry of the Li~[ cation. 




FIG. 6. Equilibrium geometry of the Lig cation. 

The lowest energy structures of LiJ and Li^ , i.e. 
the fragments corresponding to the relevant dissociation 
paths of Liii, are shown in Fig.(||) and Fig.(||), respec- 
tively. Both structures have Gi„ symmetry. By evalu- 
ating the dissociation energies for the two channels, ac- 
cording to the definitions ( |l4| ) and (|l5|), we find that the 
dissociation with the production of a monomer is ener- 
getically slightly favored over that producing a dimer. In 
fact, our calculated values for these dissociation energies 
are D+^l) = l.Q7eV and = l.lleV. 

These findings seem to support the claim, as suggested 
in Ref. , that the fragmentation process is dominated 
by non-statistical effects, which may favor dissociation 
along the less energetically favored channels. 

The sophisticated calculations of the dynamics of 
highly excited cluster which are required to explain in 
detail the experimental results of Ref. |H| are of course 
beyond the applicability of ground-state theories such 
as DFT. High level ab-initio calculations of the excited 
states incorporating large scale valence electron config- 
uration interaction are probably necessary to calculate 
reliably the optical response of the clusters and thus to 
predict the dissociation pathways. We believe, however, 
that informations about the structure and energetics such 



as those reported here may provide a useful input for 
more refined calculations where a realistic description of 
the excitation process and the way in which the energy 
is redistributed among the excited states are treated. 

V. SUMMARY AND CONCLUSIONS 

We presented a method to perform ab-initio electronic 
structure calculations in real-space by using a local and 
accurate discretization of the electron-ion Hamiltonian 
within DFT. We use optimized pseudopotentials partic- 
ularly suited for calculations in real-space. The force 
acting on the ions are calculated accurately and the fea- 
sibility of ab-initio Molecular Dynamics is demonstrated. 
At variance with a recently proposed Mehrstellen-MG 
scheme, our method implements a FMG schedule that 
uses a set of coarser grid where the KS equations are 
solved approximately. We found that, for cluster calcula- 
tions, our code has performances comparable or superior 
to the most efficient PW codes available nowadays. 

The method has been applied to ab-initio MD simu- 
lations with large time steps without loss of accuracy in 
the total energy during the simulations. As for many 
real-space schemes proposed recently, one of the major 
advantages of this FMG scheme is that it can be read- 
ily adapted to run efficiently on parallel computer ar- 
chitectures: all operations, except for the orthogonaliza- 
tion of the KS orbitals, are short-ranged. For the simple 
molecules and clusters investigated here, the largest frac- 
tion (80- 90%) of the CPU time is spent in performing 
relaxations, which execute very efficiently on computers 
with parallel architecture. We have studied with this 
method the energetics of fragmentation of Lif x clusters 
and confirmed the importance of the statistical effects in 
their dissociation. 
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